********************************************************************************
******************************Column 4******************************************
********************************************************************************

set     more off
set     matsize 800
clear   all
global  dpath = "[Userpath]\data and code"

use "$dpath\data\firmdata.dta",clear
sort countycode year

*1.County-level data

foreach var in output wastewater cod nh3n gas so2 dust {
	bysort countycode year: egen `var's = sum(`var')
	drop `var'
	rename `var's `var'
	gen l`var' = log(`var')
}
duplicates drop countycode year, force

gen llgdpC = log(lgdpC)
drop if llgdpC == .

gen llgdp = log(lgdp)
drop if llgdp == .

*2. Remove outliers
winsor2 lwastewater, replace cuts(0.5 99.5) trim
winsor2 lcod, replace cuts(0.5 99.5) trim
winsor2 lgas, replace cuts(0.5 99.5) trim
winsor2 lso2, replace cuts(0.5 99.5) trim
winsor2 ldust, replace cuts(0.5 99.5) trim

*3. Generate variables
xtset countycode year

gen manualupstream = 1 if manualyear != 33000
replace manualupstream = 0 if manualupstream == .

gen treated = 1 if autoupstream == 1 & year >= autoyear
replace treated = 0 if treated == .

gen treatedM = 1 if manualupstream == 1 & year >= manualyear
replace treatedM = 0 if treatedM == .

cd "$dpath\results"

qui tab year, g(t_)
g time=year-1997

xi i.countycode*time
qui foreach v of var _IcouXt* {
 gen name`v' = `v'
}
drop _Icou*

drop time

g post=(year>=autoyear & autoyear != 33000)
g trend=(year-autoyear)*post

xtset,clear

g prtt_lwastewater = .
g prtt_lcod = .

cap program drop pro1

program pro1
	reghdfe lwastewater name_IcouXt* t_1-t_13 treatedM llgdpC if post == 0,absorb(countycode)
	foreach m of var name_IcouXt* {
		cap replace prtt_lwastewater = _b[`m'] * trend if `m' > 0
	}
reghdfe lwastewater treated treatedM llgdpC prtt_lwastewater, absorb(countycode citycode#year) 
end

bootstrap,reps(500) seed(10101) nodots cluster(countycode):pro1
estimates store reg1

esttab reg1 using TableAIX.rtf,append b(%9.3f) se scalar(N r2_a) star(* 0.1 **  0.05  *** 0.01) title(Table AIX Column 4) keep(treated)

********************************************************************************
******************************Column 5/6****************************************
********************************************************************************
set     more off
set     matsize 800
clear   all
global  dpath = "G:\水污染与健康\data and code"

use "$dpath\data\deathratecontrol.dta",clear
cd "$dpath\results"

winsor2 lifeexpectancy, replace cuts(1 99) trim

qui tab year, g(t_)
g time=year-1990

xi i.code*time
qui foreach v of var _IcodXt* {
 gen name`v' = `v'
}
drop _Icod*
drop time

g post=(year>=autoyear & autoyear != 3300)
g trend=(year-autoyear)*post


g prtt_totaldigestivedr = .
g prtt_lifeexpectancy = .

cap program drop pro1
cap program drop pro2

xtset,clear

program pro1
	reghdfe totaldigestivedr name_IcodXt* t_1-t_19 manual lngdp lnbed if post == 0,absorb(code)
	foreach m of var name_IcodXt* {
		cap replace prtt_totaldigestivedr = _b[`m'] * trend if `m' > 0
	}
reghdfe totaldigestivedr auto upstreamauto manual lngdp lnbed prtt_totaldigestivedr, absorb(code year) 
end
bootstrap,reps(500) seed(10101) nodots cluster(code):pro1
estimates store reg1

program pro2
	reghdfe lifeexpectancy name_IcodXt* t_1-t_19 manual lngdp lnbed if post == 0,absorb(code)
	foreach m of var name_IcodXt* {
		cap replace prtt_lifeexpectancy = _b[`m'] * trend if `m' > 0
	}
reghdfe lifeexpectancy auto upstreamauto manual lngdp lnbed prtt_lifeexpectancy, absorb(code year) 
end
bootstrap,reps(500) seed(10101) nodots cluster(code):pro2
estimates store reg2

esttab reg1 reg2 using TableAIX.rtf,append b(%9.3f) se scalar(N r2_a) star(* 0.1 **  0.05  *** 0.01) title(Table AIX Column 5/6) keep(auto)
